#include<math.h>
#include"Project.cuh"
#include <iostream>
using namespace NNU::DGGS;
using namespace NNU::DGGS::Base;
using namespace NNU::DGGS::Kernel::Projection;

#pragma region paraments

//菱形块编号 三角形编号 邻接三角形序号 keep X平移量 Y平移量 旋转系数
__constant__  double _vertTable[504] = {
        // vert 0
        2, 1, 0, 1, -M_HALF, -M_SIN60, 3,
        1, 0, 1, 1, -M_ONE, M_ZERO, 2,
        5, 4, 2, 1, -M_HALF, M_SIN60, 1,
        -1, -1, 3, 0, -M_HALF, M_SIN60, 1,
        4, 3, 4, 1, M_ONE, M_ZERO, -1,
        3, 2, 5, 1, M_HALF, -M_SIN60, -2,

        // vert 1
        1, 0, 0, 1, M_ZERO, M_ZERO, 1,
        1, 5, 1, 1, -M_HALF, -M_SIN60, 4,
        10, 14, 2, 1, -M_HALF, M_SIN60, 1,
        -1, -1, 3, 0, -M_HALF, M_SIN60, 1,
        5, 9, 4, 1, M_ZERO, M_ZERO, 3,
        5, 4, 5, 1, M_ONE, M_ZERO, 0,

        // vert 2
        2, 1, 0, 1, M_ZERO, M_ZERO, 1,
        2, 6, 1, 1, -M_HALF, -M_SIN60, 4,
        10, 10, 2, 1, -M_HALF, M_SIN60, 1,
        -1, -1, 3, 0, -M_HALF, M_SIN60, 1,
        1, 5, 4, 1, M_ZERO, M_ZERO, 3,
        1, 0, 5, 1, M_ONE, M_ZERO, 0,

        // vert 3
        3, 2, 0, 1, M_ZERO, M_ZERO, 1,
        3, 7, 1, 1, -M_HALF, -M_SIN60, 4,
        7, 11, 2, 1, -M_HALF, M_SIN60, 1,
        -1, -1, 3, 0, -M_HALF, M_SIN60, 1,
        2, 6, 4, 1, M_ZERO, M_ZERO, 3,
        2, 1, 5, 1, M_ONE, M_ZERO, 0,

        // vert 4
        4, 3, 0, 1, M_ZERO, M_ZERO, 1,
        4, 8, 1, 1, -M_HALF, -M_SIN60, 4,
        8, 12, 2, 1, -M_HALF, M_SIN60, 1,
        -1, -1, 3, 0, -M_HALF, M_SIN60, 1,
        3, 7, 4, 1, M_ZERO, M_ZERO, 3,
        3, 2, 5, 1, M_ONE, M_ZERO, 0,

        // vert 5
        5, 4, 0, 1, M_ZERO, M_ZERO, 1,
        5, 9, 1, 1, -M_HALF, -M_SIN60, 4,
        9, 13, 2, 1, -M_HALF, M_SIN60, 1,
        -1, -1, 3, 0, -M_HALF, M_SIN60, 1,
        4, 8, 4, 1, M_ZERO, M_ZERO, 3,
        4, 3, 5, 1, M_ONE, M_ZERO, 0,

        // vert 6
        6, 10, 0, 1, M_ZERO, M_ZERO, 1,
        6, 15, 1, 1, -M_HALF, -M_SIN60, 4,
        10, 19, 2, 1, M_ZERO, M_ZERO, -1,
        10, 14, 3, 1, -M_HALF, M_SIN60, 2,
        -1, -1, 4, 0, -M_HALF, M_SIN60, 1,
        1, 5, 5, 1, M_HALF, -M_SIN60, 4,

        // vert 7
        7, 11, 0, 1, M_ZERO, M_ZERO, 1,
        7, 16, 1, 1, -M_HALF, -M_SIN60, 4,
        6, 15, 2, 1, M_ZERO, M_ZERO, -1,
        6, 10, 3, 1, -M_HALF, M_SIN60, 2,
        -1, -1, 4, 0, -M_HALF, M_SIN60, 1,
        2, 6, 5, 1, M_HALF, -M_SIN60, 4,

        // vert 8
        8, 12, 0, 1, M_ZERO, M_ZERO, 1,
        8, 17, 1, 1, -M_HALF, -M_SIN60, 4,
        7, 16, 2, 1, M_ZERO, M_ZERO, -1,
        7, 11, 3, 1, -M_HALF, M_SIN60, 2,
        -1, -1, 4, 0, -M_HALF, M_SIN60, 1,
        3, 7, 5, 1, M_HALF, -M_SIN60, 4,

        // vert 9
        9, 13, 0, 1, M_ZERO, M_ZERO, 1,
        9, 18, 1, 1, -M_HALF, -M_SIN60, 4,
        8, 17, 2, 1, M_ZERO, M_ZERO, -1,
        8, 12, 3, 1, -M_HALF, M_SIN60, 2,
        -1, -1, 4, 0, -M_HALF, M_SIN60, 1,
        4, 8, 5, 1, M_HALF, -M_SIN60, 4,

        // vert 10
        10, 14, 0, 1, M_ZERO, M_ZERO, 1,
        10, 19, 1, 1, -M_HALF, -M_SIN60, 4,
        9, 18, 2, 1, M_ZERO, M_ZERO, -1,
        9, 13, 3, 1, -M_HALF, M_SIN60, 2,
        -1, -1, 4, 0, -M_HALF, M_SIN60, 1,
        5, 9, 5, 1, M_HALF, -M_SIN60, 4,

        // vert 11
        8, 17, 0, 1, -M_HALF, -M_SIN60, 3,
        9, 18, 1, 1, -M_ONE, M_ZERO, 2,
        10, 19, 2, 1, -M_HALF, M_SIN60, 1,
        6, 15, 3, 1, M_HALF, M_SIN60, 0,
        -1, -1, 4, 0, M_ZERO, M_ZERO, 0,
        7, 16, 5, 1, M_HALF, -M_SIN60, -2
};


//lon lat 方位角 旋转系数 平移量X 平移量Y 菱形块ID
__constant__ double _triParaments[140] = {
        -2.5132741229795865, 0.91843818699424362, 0.000000, 1, 0.000000, 0.000000, 1,
        -1.2566370614897933, 0.91843818696277224, 0.000000, 1, 0.000000, 0.000000, 2,
        -0.0000000000000000, 0.91843818699424362, 0.000000, 1, 0.000000, 0.000000, 3,
        1.2566370614897928, 0.91843818696277224, -0.000000, 1, 0.000000, 0.000000, 4,
        2.5132741228, 0.91843818693130097, -0.000000, 1, 0.000000, 0.000000, 5,
        -2.5132741229452886, 0.18871053073216107, 3.1415926535, 4, -0.500000, -0.866025, 1,
        -1.2566370614897933, 0.18871053072473168, 3.1415926535, 4, -0.500000, -0.866025, 2,
        -0.0000000000000000, 0.18871053073216107, 3.1415926535, 4, -0.500000, -0.866025, 3,
        1.256637061489793, 0.18871053072473171, 3.1415926535, 4, -0.500000, -0.866025, 4,
        2.5132741228342979, 0.18871053071730234, 3.1415926535, 4, -0.500000, -0.866025, 5,
        -1.8849555921897931, -0.18871053082595737, -0.000000, 1, 0.000000, 0.000000, 6,
        -0.62831853076204569, -0.1887105308296721, -0.000000, 1, 0.000000, 0.000000, 7,
        0.62831853076204558, -0.1887105308296721, 0.000000, 1, 0.000000, 0.000000, 8,
        1.8849555921897931, -0.18871053082595743, -0.000000, 1, 0.000000, 0.000000, 9,
        3.1415926535211973, -0.18871053082595737, -0.000000, 1, 0.000000, 0.000000, 10,
        -1.8849555922192975, -0.91843818702070867, -3.1415926535, 4, -0.500000, -0.866025, 6,
        -0.62831853076313138, -0.91843818705317659, -3.1415926535, 4, -0.500000, -0.866025, 7,
        0.62831853076313138, -0.91843818705317659, 3.1415926535, 4, -0.500000, -0.866025, 8,
        1.8849555922192978, -0.91843818702070867, 3.1415926535, 4, -0.500000, -0.866025, 9,
        -3.1415926535897931, -0.91843818701036806, -3.1415926535, 4, -0.500000, -0.866025, 10
};

#pragma endregion

/// <summary>
/// 初始化
/// </summary>
__DGGS_CPGPU__ void Project::init()
{

}



/********************************************************格网编码转经纬度坐标*******************************************************/
#pragma region RC2LBCenter
/// <summary>
/// 行列编码转经纬度
/// </summary>
__DGGS_CPGPU__ void Project::RC2LBCenter(DGGSGridType topo, int basePartition, long long int row, long long int col, int level, double &lon, double &lat) {
    double result[2] = {0.0, 0.0};

    double cent[2] = {0.0, 0.0};
    invQuantify(topo,row, col, cent);

    double trans[2] = {0.0, 0.0};
    double rotDeg = 0.0;
    double scaleFac = pow(2, level);

    if (topo == DGGSGridType::DGGS_ISEA4H)
    {
        trans[0]=0.0,trans[1]= 0.0;
        scaleFac = (double)(pow(2, level));
    }
    else if (topo == DGGSGridType::DGGS_ISEA4T)
    {
        trans[0]=-sqrtThree / 2,trans[1]= -0.5l;
        scaleFac = (double)(pow(2, level) * sqrtThree);
    }
    else if (topo == DGGSGridType::DGGS_ISEA4D)
    {
        trans[0]=-0.5l / 2,trans[1]= -offsetY;
        scaleFac = pow(2, level);
    }
    else
    {
        return;
    }


    convertTypeAddress(cent, scaleFac, trans, rotDeg, basePartition, result);
    lon = result[0] ;
    lat = result[1] ;
}

/// <summary>
/// 计算两点在投影平面上的距离
/// </summary>
__DGGS_CPGPU__ double Project::getSnyderDist(Pt2D &g1, const Pt2D &g2, bool rads)
{
    double bigC = fabs(g2.x - g1.x);

    if (bigC > M_PI) {
        double lon1 = g1.x;
        if (lon1 < 0.0L) lon1 += 2.0L * M_PI;
        double lon2 = g2.x;
        if (lon2 < 0.0L) lon2 += 2.0L * M_PI;

        bigC = fabs(lon2 - lon1);
    }

    double b = M_PI / 2 - g1.y;
    double a = M_PI / 2 - g2.y;


    double cosc = cos(a) * cos(b) + sin(a) * sin(b) * cos(bigC);

    if (cosc > 1.0L) cosc = 1.0L;
    if (cosc < -1.0L) cosc = -1.0L;

    double retVal = acos(cosc);

    if (!rads) retVal *= radToDeg;

    return retVal;

}

/// <summary>
/// 计算格子落在顶点坐标系的哪个三角形上
/// </summary>
__DGGS_CPGPU__ int Project::calSubtriangle(const double &x, const double &y) {

    double xs = sqrtThree * x;

    double xpp = xs + tol,
            xmp = -xs + tol,
            xpm = xs - tol,
            xmm = -xs - tol;

    if (y >= xmm && y > xpp) {
        return 0;
    }

    if ((fabs(y) <= tol && fabs(x) <= tol) ||
        (y <= xpp && y >= (0.0 - tol))) {
        return 1;
    }

    if (y < (0.0 - tol) && y > xmp) {
        return 2;
    }

    if (y <= xmp && y < xpm) {
        return 3;
    }

    if (y >= xpm && y < (0.0 - tol)) {
        return 4;
    }

    if (y >= (0.0 - tol) && y < xmm) {
        return 5;
    }

    return -1;
}

/// <summary>
/// 行列号转格网坐标
/// </summary>
__DGGS_CPGPU__ void Project::invQuantify(DGGSGridType topo, long long int i, long long int j, double *add)
{
    if (topo == DGGSGridType::DGGS_ISEA4H || topo == DGGSGridType::DGGS_ISEA3H)
    {
        add[0] = (double)i - 0.5 * (double)j;
        add[1] = (double)j * 1.5 * sqrtThree / 3;
    }
    else if (topo == DGGSGridType::DGGS_ISEA4T)
    {
        add[0] = (double)i * sqrtThree - ((double)((j + 1) / 2)* sqrtThree / 2.0);
        add[1] = 1.5 * (double)(j / 2) + 0.5 * (j % 2);
    }
    else
    {
        add[0] = (double)i;
        add[1] = (double)j * offsetYDouble;

        add[0] = add[0] - skewFactor * add[1];
    }
}

/// <summary>
/// 施耐德反算
/// </summary>
__DGGS_CPGPU__ void Project::snyderInv(double *xyCoord, int triangle, double *lonlatCoord) {
    double ddazh, ph, fh, azh, azh1, dazh, h, fazh, flazh, dz, z;
    double sinLat, sinLon;
    double azh0;

    //cent
    double triLon = _triParaments[triangle * 7 + 0];
    double triLat = _triParaments[triangle * 7 + 1];

    //pt
    double ptX = xyCoord[0] * icosaEdgeLength - originXOffset;
    double ptY = xyCoord[1] * icosaEdgeLength - originYOffset;

    ddazh = _triParaments[triangle * 7 + 2];


    if ((fabs(ptX) < PRECISION) && (fabs(ptY) < PRECISION)) {
        lonlatCoord[0] = triLon;
        lonlatCoord[1] = triLat;
    }
    else
    {
        ph = sqrt(ptX * ptX + ptY * ptY);
        azh1 = atan2(ptX, ptY);

        if (azh1 < 0.0L)
            azh1 = azh1 + 2 * M_PI;
        azh0 = azh1;
        if ((azh1 > 120.0L * degToRad) && (azh1 <= 240.0L * degToRad))
            azh1 = azh1 - 120.0L * degToRad;
        if (azh1 > 240.0L * degToRad)
            azh1 = azh1 - 240.0L * degToRad;
        azh = azh1;

        if (fabs(azh1) > PRECISION)
        {
             double agh = r1Squared * tanDh * tanDh / (2.0L * (1.0L / tan(azh1) + cot30));
            dazh = 1.0;

            while (fabs(dazh) > PRECISION) {

                h = acos(sin(azh) * sinGh * cosDh - cos(azh) * cosGh);
                fazh = agh - azh - GH - h + M_PI;
                flazh = ((cos(azh) * sinGh * cosDh + sin(azh) * cosGh) / sin(h)) - 1.0;
                dazh = -fazh / flazh;
                azh = azh + dazh;
            }
        }
        else
            azh = azh1 = 0.0;


        dz = atan2(tanDh, cos(azh) + cot30 * sin(azh));
        fh = tanDh / (2.0 * (cos(azh1) + cot30 * sin(azh1)) * sin(dz / 2.0));
        z = 2.0 * asin(ph / (2.0 * r1Value * fh));
        if ((azh0 >= 120 * degToRad) && (azh0 < 240.0 * degToRad)) azh = azh + 120 * degToRad;
        if (azh0 >= 240.0 * degToRad) azh = azh + 240.0 * degToRad;

        azh += ddazh;

        while (azh <= -M_PI) azh += M_2PI;
        while (azh > M_PI) azh -= M_2PI;

        sinLat = sin(triLat) * cos(z) + cos(triLat) * sin(z) * cos(azh);
        if (sinLat > M_ONE) sinLat = M_ONE;
        if (sinLat < -M_ONE) sinLat = -M_ONE;
        lonlatCoord[1] = asin(sinLat);

        if (fabs(fabs(lonlatCoord[1]) - PIHalf) < PRECISION)
        {
            lonlatCoord[1] = (lonlatCoord[1] > M_ZERO) ? PIHalf : -PIHalf;
            lonlatCoord[0] = M_ZERO;
        }
        else
        {
            sinLon = sin(azh) * sin(z) / cos(lonlatCoord[1]);
            double cosLon = (cos(z) - sin(triLat) * sin(lonlatCoord[1])) /
                                 cos(triLat) / cos(lonlatCoord[1]);
            if (sinLon > M_ONE) sinLon = M_ONE;
            if (sinLon < -M_ONE) sinLon = -M_ONE;
            if (cosLon > M_ONE) cosLon = M_ONE;
            if (cosLon < -M_ONE) cosLon = -M_ONE;
            lonlatCoord[0] = triLon + asin(sinLon);
            lonlatCoord[0] = triLon + atan2(sinLon, cosLon);
            if (lonlatCoord[0] <= -M_PI) lonlatCoord[0] += M_2PI;
            if (lonlatCoord[0] >= M_PI) lonlatCoord[0] -= M_2PI;
        }
    }

    lonlatCoord[0]*= radToDeg;
    lonlatCoord[1]*= radToDeg;
}

/// <summary>
/// 编码转经纬度
/// </summary>
__DGGS_CPGPU__ void Project::convertTypeAddress(double *cent, double scaleFac, double *trans, double rotDeg, int quadNum,
                              double *result) {

    double tempRes[2] = {0.0, 0.0};

    double x_temp = cent[0]-trans[0];
    double y_temp = cent[1]-trans[1];
    x_temp /= scaleFac;
    y_temp /= scaleFac;

    tempRes[0] = x_temp;
    tempRes[1] = y_temp;

    result[0] = x_temp;
    result[1] = y_temp;

    int subTri = calSubtriangle(x_temp, y_temp);


    if (_vertTable[quadNum * 42 + subTri * 7 + 3] == 1.0)
    {
        double rotateAng = _vertTable[quadNum * 42 + subTri * 7 + 6] * -60.0;
        while (rotateAng < M_ZERO) rotateAng += 360.0L;
        while (rotateAng >= 360.0L) rotateAng -= 360.0L;
        rotateAng *= degToRad;

        x_temp += _vertTable[quadNum * 42 + subTri * 7 + 4];
        y_temp += _vertTable[quadNum * 42 + subTri * 7 + 5];

        tempRes[0] = x_temp * cos(rotateAng) - y_temp * sin(rotateAng);
        tempRes[1] = x_temp * sin(rotateAng) + y_temp * cos(rotateAng);

        int triangle = (int)(_vertTable[quadNum * 42 + subTri * 7 + 1]);

        double lonlatCoord[2] = {0.0, 0.0};

        //这一步加上算不出来结果
        snyderInv(tempRes, triangle, lonlatCoord);

        if (fabs(lonlatCoord[1]) >= (90.0 - PRECISION)) {
            lonlatCoord[0] = 0.0;
        }
        else {
            while (lonlatCoord[0] <= -180.0) {
                lonlatCoord[0] = (lonlatCoord[0]+ 360.0) ;
            }
            while (lonlatCoord[0]> 180.0) {
                lonlatCoord[0] = (lonlatCoord[0] - 360.0);
            }
        }

        result[0] = lonlatCoord[0];
        result[1] = lonlatCoord[1];

    }
}

#pragma endregion
/********************************************************经纬度坐标转格网编码*******************************************************/
#pragma region LB2RC
/// <summary>
/// 经纬度转行列坐标
/// </summary>
__DGGS_CPGPU__ bool Project::LB2RC(double lon, double lat, int level, DGGSGridType topo, int &basePartition, long long int &row, long long int &col)
{
    //经纬度转换为弧度制
    Pt2D coord;
    coord.y = lat * degToRad;
    coord.x = lon * degToRad;

    //从球体坐标转换为平面坐标
    SurfacePt gridpt = this->snyderFwd(coord);

    //缩放系数
    double scaleFac = 0;
    //平移量
    Pt2D trans;
    trans.x=0;
    trans.y=0;
    if (topo == DGGSGridType::DGGS_ISEA4H)
    {
        scaleFac = pow(2, level);
    }
    else if (topo == DGGSGridType::DGGS_ISEA4T)
    {
        scaleFac =  (double)(pow(2, level) * sqrtThree);
        trans.x=-sqrtThreeHalf;
        trans.y=-1.0l / 2;

    }
    else if (topo == DGGSGridType::DGGS_ISEA4D)
    {
        scaleFac = pow(2, level);
        trans.x=-1.0l / 4;
        trans.y=-sqrtThreeHalf / 2;
    }
    else
    {
        return false;
    }

    //对平面坐标进行旋转、平移和缩放，得到格网行列号
    RCPt2D newCoord;
    basePartition = planeCoords2RC(gridpt, scaleFac, trans,topo, newCoord);
    row = (newCoord.i);
    col = (newCoord.j);

    /*
     * 修复六边形投影误差错误
     * */
    if (topo == DGGSGridType::DGGS_ISEA4H)
    {
        Pt2D tmpCoord;
        tmpCoord.x = gridpt.x;
        tmpCoord.y = gridpt.y;

        rotate(_triParaments[gridpt.triangle * 7 + 3] * 60.0L, tmpCoord);

        tmpCoord.x -= _triParaments[gridpt.triangle * 7 + 4];
        tmpCoord.y -= _triParaments[gridpt.triangle * 7 + 5];

        tmpCoord = scale(tmpCoord, scaleFac);
        tmpCoord.x += trans.x;
        tmpCoord.y += trans.y;

        //repair(topo, tmpCoord.x, tmpCoord.y,level,row,col);
    }

    if(col<0&&row<0)
    {
        col=0;
        row=0;
    }
    else {
        if (row < 0) {
            if (basePartition < 6) {
                basePartition = basePartition - 1;
                if (basePartition < 1)
                    basePartition = 5;
                row = (long long int)(pow(2, level) - 1 - col);
                col = (long long int)(pow(2, level) - 1);
            } else {
                basePartition = basePartition - 6;
                if (basePartition < 1)
                    basePartition = 10;
                col = col;
                row =(long long int)(pow(2, level) - 1);
            }
        }
        else if (row >= pow(2, level)) {
            if (basePartition < 6) {
                basePartition = basePartition + 5;
                row = 0;
                col = col;
            } else {
                basePartition = basePartition + 1;
                if (basePartition > 10)
                    basePartition = 6;
                row = (long long int)(pow(2, level) - col);
                col = 0;
                if (row >= pow(2, level)) {
                    basePartition = 11;
                    row = 0;
                    col = 0;
                }
            }
        }
        if (topo == DGGSGridType::DGGS_ISEA4T)
            level++;
        if (col < 0) {
            if (basePartition < 6) {
                basePartition = basePartition + 4;
                if (basePartition < 6)
                    basePartition = 10;
                row = row;
                col = (long long int)pow(2, level) - 1;
            } else {
                basePartition = basePartition + 1;
                if (basePartition > 10)
                    basePartition = 6;
                col = (long long int)pow(2, level) - 1 - col;
                row = (long long int)pow(2, level) - 1;
            }
        }
        else if (col >= pow(2, level)) {
            if (basePartition < 6) {
                basePartition = basePartition + 1;
                if (basePartition > 5)
                    basePartition = 1;
                col = (long long int)pow(2, level) - row;
                row = 0;
                if (col >= pow(2, level)) {
                    basePartition = 0;
                    row = 0;
                    col = 0;
                }
            } else {
                basePartition = basePartition - 4;
                if (basePartition > 5)
                    basePartition = 1;
                row = row;
                col = 0;
            }
        }
    }
    //将行列号转换为格元编号,南北极点情况特殊，需要分类处理
    return true;
}


/// <summary>
/// 球面坐标转平面坐标
/// </summary>
__DGGS_CPGPU__ Pt2D Project::sphere2Plane(const Pt2D &geoVect, int nTri) {
     double ph, fh, azh, azh1, dazh, h, dz, z, azh0, ag, cosAzh, sinAzh;
    Pt2D Planevect;

    //中心点
     double tri_lon = _triParaments[nTri * 7 + 0];
     double tri_lat = _triParaments[nTri * 7 + 1];

     double cosLat = cos(geoVect.y);
     double sinLat = sin(geoVect.y);

    dazh = _triParaments[nTri * 7 + 2];

     double tmp = sin(tri_lat) * sinLat + cos(tri_lat) * cosLat *
                                              cos(geoVect.x - tri_lon);

    if (tmp > 1.0) tmp = 1.0;
    if (tmp < -1.0) tmp = -1.0;
    z = acos(tmp);

    azh = atan2(cosLat * sin(geoVect.x - tri_lon),
                cos(tri_lat) * sinLat - sin(tri_lat) * cosLat * cos(geoVect.x - tri_lon)) - dazh;

    if (azh < 0.0) azh = azh + 2.0 * M_PI;
    azh0 = azh;
    if ((azh >= 120.0 * degToRad) && (azh <= 240.0 * degToRad)) azh -= 120.0 * degToRad;
    if (azh > 240.0 * degToRad) azh -= 240.0 * degToRad;

    cosAzh = cos(azh);
    sinAzh = sin(azh);

    dz = atan2(tanDh, cosAzh + cot30 * sinAzh);

    h = acos(sinAzh * sinGh * cosDh - cosAzh * cosGh);
    ag = azh + GH + h - 180.0 * degToRad;
    azh1 = atan2(2.0 * ag, r1Squared * tanDh * tanDh - 2.0 * ag * cot30);
    fh = tanDh / (2.0 * (cos(azh1) + cot30 * sin(azh1)) * sin(dz/2.0));
    ph = 2.0 * r1Value * fh * sin(z/2.0);

    if ((azh0 >= 120.0 * degToRad) && (azh0 < 240.0 * degToRad)) azh1 += 120.0* degToRad;
    if (azh0 >= 240.0 * degToRad) azh1 += 240 * degToRad;

    Planevect.x = (ph * sin(azh1) + originXOffset) / icosaEdgeLength;
    Planevect.y = (ph * cos(azh1) + originYOffset) / icosaEdgeLength;
    return Planevect;

}

/// <summary>
/// 计算点属于哪个三角面
/// </summary>
__DGGS_CPGPU__ int Project::computeSubTriangle(const Pt2D &pt) {
    // 从第一个面开始计算，逐个比较
    int minFace = 0;

    Pt2D coord;
    coord.x= _triParaments[0];
    coord.y= _triParaments[1];


    double minDist = getSnyderDist(coord, pt, true);

    for (int i = 1; i < 20; i++) {
        Pt2D triCen;
        triCen.x= _triParaments[7 * i];
        triCen.y= _triParaments[7 * i + 1];
        double dist = getSnyderDist(triCen, pt, true);
        if (dist < minDist) {
            minDist = dist;
            minFace = i;
        }
    }
    return minFace;
}

/// <summary>
/// 施奈德正算
/// </summary>
__DGGS_CPGPU__ SurfacePt Project::snyderFwd(const Pt2D &ll) {
    SurfacePt gridpt;

    gridpt.triangle = computeSubTriangle(ll);

    if (gridpt.triangle < 0) {
        gridpt.x = 0.0l;
        gridpt.y = 0.0l;

        return gridpt;
    }

    Pt2D pt = sphere2Plane(ll, gridpt.triangle);
    gridpt.x = pt.x;
    gridpt.y = pt.y;
    return gridpt;

}

/// <summary>
/// 旋转算法
/// </summary>
__DGGS_CPGPU__ void Project::rotate(double degrees, Pt2D &point) {
    while (degrees < 0) degrees += 360.0L;
    while (degrees >= 360.0L) degrees -= 360.0L;

    double rotAng = degrees * degToRad;
    double cosAng = cos(rotAng);
    double sinAng = sin(rotAng);

    double x = point.x;
    double y = point.y;
    point.x = x * cosAng - y * sinAng;
    point.y = x * sinAng + y * cosAng;

}

/// <summary>
/// 缩放算法
/// </summary>
__DGGS_CPGPU__ Pt2D Project::scale(const Pt2D addIn, double scaleFac) {
    Pt2D newAdd(addIn);
    newAdd.x *= scaleFac;
    newAdd.y *= scaleFac;
    return newAdd;
}


/// <summary>
/// 经纬度转行列坐标
/// </summary>
__DGGS_CPGPU__ int Project::planeCoords2RC(const SurfacePt &gridpt, double scaleFac, Pt2D trans, DGGSGridType topo, RCPt2D &result) {
    int quadNum = -1;

    Pt2D tmpCoord;
    tmpCoord.x = gridpt.x;
    tmpCoord.y = gridpt.y;

    rotate(_triParaments[gridpt.triangle * 7 + 3] * 60.0L, tmpCoord);

    tmpCoord.x -= _triParaments[gridpt.triangle * 7 + 4];
    tmpCoord.y -= _triParaments[gridpt.triangle * 7 + 5];

    quadNum = (int)(_triParaments[gridpt.triangle * 7 + 6]);

    tmpCoord = scale(tmpCoord, scaleFac);
    tmpCoord.x += trans.x;
    tmpCoord.y += trans.y;

    if (topo == DGGSGridType::DGGS_ISEA4H)
    {
        result = HexQuantify(tmpCoord);
    }
    else if (topo == DGGSGridType::DGGS_ISEA4D)
    {
        result = rhoQuantify(tmpCoord);
    }
    else if (topo == DGGSGridType::DGGS_ISEA4T)
    {
        result = triQuantify(tmpCoord);
    }
    return quadNum;

}

/// <summary>
/// 菱形格子量化坐标
/// </summary>
__DGGS_CPGPU__ RCPt2D Project::rhoQuantify(const Pt2D &point) {
    Pt2D tmp;
    tmp.x = point.x;
    tmp.y = point.y;
    tmp.x = tmp.x + skewFactor * tmp.y;
    double x = offsetYDouble;

    RCPt2D result{(long long int)tmp.x,
                  (long long int)(tmp.y / x)};
    return result;
}

/// <summary>
/// 三角形格子量化坐标
/// </summary>
__DGGS_CPGPU__ RCPt2D Project::triQuantify(const Pt2D &point)
{
    Pt2D p=point;
    p.x += sqrtThree / 2.0L;
    p.y += 0.5L;

    Pt2D pp = {static_cast<double>(p.x + p.y / sqrtThree), p.y};

    long long int i = static_cast<long long int>(floor(pp.x / sqrtThree));
    long long int jTmp = static_cast<long long int>(floor(2.0 / 3.0 * pp.y));

    Pt2D ppp = {static_cast<double>(p.x - p.y / sqrtThree), p.y};
    long long int k = static_cast<long long int>(floor(ppp.x / sqrtThree));

    long long int j = jTmp * 2;
    if (jTmp % 2) {
        if (!((i + j + k) % 2)) j++;
    } else {
        if ((i + j + k) % 2) j++;
    }

    RCPt2D result = {i, j};
    return result;

}

/// <summary>
/// 六边形格子量化坐标
/// </summary>
__DGGS_CPGPU__ RCPt2D Project::HexQuantify(const Pt2D &point)
{
    double a1, a2;
    double x1, x2;
    long long int m1, m2;
    double r1, r2;
    RCPt2D add;

    a1 = point.x;
    a2 = point.y;


    x2 = a2 / sqrtThreeHalf;
    x1 = a1 + x2 / 2.0L;


    m1 = static_cast<long long int>(x1);
    m2 = static_cast<long long int>(x2);


    r1 = x1 - m1;
    r2 = x2 - m2;

    if(r1<0&&r2<0)
    {
        r1=1.0+r1;
        r2=1.0+r2;
        m1=m1-1;
        m2=m2-1;
    }
    else if(r1<0&&r2>0)
    {
        r1=1.0+r1;
        m1=m1-1;
    }
    else if(r1>0&&r2<0)
    {
        r2=1.0+r2;
        m2=m2-1;
    }

    if(r1<1.0/3.0)
    {
        if(r2<(r1+1.0)/2.0)
        {
            add.i=m1;
            add.j=m2;
        }
        else
        {
            add.i=m1;
            add.j=m2+1;
        }
    }
    else if(r1<1.0/2.0)
    {
        if(r2>2.0*r1)
        {
            add.i=m1;
            add.j=m2+1;
        }
        else if(r2<1.0-r1)
        {
            add.i=m1;
            add.j=m2;
        }
        else
        {
            add.i=m1+1;
            add.j=m2+1;
        }

    }
    else if(r1<2.0/3.0) {
        if (r2 > 1.0 - r1)
        {
            add.i=m1+1;
            add.j=m2+1;
        }
        else if(r2<2*r1-1)
        {
            add.i=m1+1;
            add.j=m2;
        }
        else
        {
            add.i=m1;
            add.j=m2;
        }
    }
    else
    {
        if(r2>r1/2.0)
        {
            add.i=m1+1;
            add.j=m2+1;
        }
        else
        {
            add.i=m1+1;
            add.j=m2;
        }
    }

    return add;
}

#pragma endregion
/********************************************************经纬度坐标转格网编码*******************************************************/
#pragma region RC2LBVertex

__DGGS_CPGPU__ bool Project::RC2LBVertex(int basePartition, long long int row, long long int col, int level, DGGSGridType topo, double* result)
{
    //格元中心点
    double cent[2] = {0.0, 0.0};
    invQuantify(topo, row, col, cent);
    //平移、旋转、缩放参数
    double trans[2]= {0.0, 0.0};;
    double rotDeg;
    double scaleFac;

    double temp[2];

    if (topo == DGGSGridType::DGGS_ISEA4H)
    {
        GeoCoord verts[6];

        if (row == 0)
            invQuantify(topo, row, (col + 1) / 2, cent);
        else
            invQuantify(topo, row, (col / 2), cent);

        double r = sqrtThree / 3;
        double r2 = r / 2;
        trans[0] = 0, trans[1] = 0;
        rotDeg = 0.0L;
        scaleFac = pow(2, level);
        if (basePartition == 0)
        {
            invQuantify(topo, 0, (long long int) scaleFac - 1, cent);
            double currentCent[2];
            currentCent[0] = cent[0] - 0.5;
            currentCent[1] = cent[1] + r2;
            if (col >= 3)
                convertTypeAddress(currentCent, scaleFac, trans, rotDeg, (int)col - 2, temp);
            else
                convertTypeAddress(currentCent, scaleFac, trans, rotDeg, (int)col + 3, temp);
        }
        else if (basePartition == 11)
        {
            invQuantify(topo, (long long int)scaleFac - 1, 0, cent);
            double currentCent[2];
            currentCent[0] = cent[0] + 0.5;
            currentCent[1] = cent[1] - r2;
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, 10 - (int)col, temp);
        }
        else
        {
            double currentCent[2];
            currentCent[0] = cent[0];
            currentCent[1] = cent[1] + r;
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, basePartition, temp);
        }
        verts[0].lon=temp[0];
        verts[0].lat=temp[1];

        double xArr[5]={- 0.5,- 0.5,0,0.5,0.5};
        double yArr[5]={r2 ,- r2,- r,- r2, r2};

        for(int i=0;i<5;i++)
        {
            double currentCent[2];
            currentCent[0] = cent[0] + xArr[i];
            currentCent[1] = cent[1] + yArr[i];
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, basePartition, temp);
            verts[i+1].lon=temp[0];
            verts[i+1].lat=temp[1];
        }
        if (basePartition == 0 || basePartition == 11)
        {
            result[0] =  (double)verts[0].lon;
            result[1] =  (double)verts[0].lat;
        }
        else if (col == 0 && row == 0)
        {
            result[0] =  (double)verts[2].lon;
            result[1] =  (double)verts[2].lat;
        }
        else if (row == 0)
        {
            if (col % 2 == 1)
            {
                result[0] =  (double)verts[3].lon;
                result[1] =  (double)verts[3].lat;
            }
            else
            {
                result[0] =  (double)verts[2].lon;
                result[1] =  (double)verts[2].lat;
            }
        }
        else
        {
            if (col % 2 != 0)
            {
                result[0] =  (double)verts[2].lon;
                result[1] =  (double)verts[2].lat;
            }
            else
            {
                result[0] =  (double)verts[3].lon;
                result[1] =  (double)verts[3].lat;
            }
        }
    }
    else if(topo == DGGSGridType::DGGS_ISEA4D)
    {
        GeoCoord vert;
        trans[0] = -0.5l / 2,trans[1] = -offsetY;
        rotDeg = 0.0;
        scaleFac = pow(2, level);
        if (basePartition == 0)
        {
            invQuantify(topo, 0, (long long int) scaleFac - 1, cent);
            double currentCent[2];
            currentCent[0] = cent[0] - rhoxOffsetComplement;
            currentCent[1] = cent[1] + rhoyOffset;
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, 2, temp);
            vert.lon = temp[0];
            vert.lat = temp[1];
        }
        else if (basePartition == 11)
        {
            invQuantify(topo, (long long int)scaleFac - 1, 0, cent);
            double currentCent[2];
            currentCent[0] = cent[0] + rhoxOffsetComplement;
            currentCent[1] = cent[1] - rhoyOffset;
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, 7, temp);
            vert.lon = temp[0];
            vert.lat = temp[1];
        }
        else
        {
            double currentCent[2];
            currentCent[0] = cent[0] - rhoxOffset;
            currentCent[1] = cent[1] - rhoyOffset;
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, basePartition, temp);
            vert.lon = temp[0];
            vert.lat = temp[1];
        }
        result[0] = vert.lon;
        result[1] = vert.lat;
    }
    else if (topo == DGGSGridType::DGGS_ISEA4T)
    {
        GeoCoord vert;
        trans[0] = -sqrtThree / 2,trans[1] = -0.5l;
        rotDeg = 0.0;
        scaleFac =  (double)(pow(2, level) * sqrtThree);

        if (basePartition == 0)
        {
            vert.lon = 0.0;
            vert.lat = 90.0;
        }
        else if (basePartition == 11)
        {
            vert.lon = 0.0;
            vert.lat = -90.0;
        }
        else
        {
            invQuantify(topo,row, col * 2 + 1,cent);
            double currentCent[2];
            currentCent[0] = cent[0];
            currentCent[1] = cent[1] - 1.0;
            convertTypeAddress(currentCent, scaleFac, trans, rotDeg, basePartition, temp);
            vert.lon = temp[0];
            vert.lat = temp[1];
        }
        result[0] = vert.lon;
        result[1] = vert.lat;
    }
    else
    {
        return false;
    }
    return true;
}

#pragma endregion